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Abstract: We previously developed the direct interaction approximation (DIA) method to 
estimate the protein-ligand binding free energy (AG). The DIA method estimates the AG 
value based on the direct van der Waals and electrostatic interaction energies between the 
protein and the ligand. In the current study, the effect of the entropy of the ligand was 
introduced with protein dynamic properties by molecular dynamics simulations, and the 
interaction between each residue of the protein and the ligand was also weighted considering 
the hydration of each residue. The molecular dynamics simulation of the apo target protein 
gave the hydration effect of each residue, under the assumption that the residues, which 
strongly bind the water molecules, are important in the protein-ligand binding. These two 
effects improved the reliability of the DIA method. In fact, the parameters used in the 
DIA became independent of the target protein. The averaged error of AG estimation was 
1.3 kcal/mol and the correlation coefficient between the experimental AG value and the 
calculated AG value was 0.75. 

Keywords: protein-ligand docking; molecular dynamics simulation; protein-ligand binding 
free energy 
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1. Introduction 

In the pharmaceutical sciences, the protein- ligand binding free energy (AG) is one of the most 
important properties of a drug compound. Despite the development of numerous docking programs 
and scoring functions to estimate the AG [1—7], the typical accuracy of AG estimation remains about 
2-3 kcal/mol [6-10]. Usually, docking scores are proportional to AG values. This low accuracy of the AG 
or docking score has contributed to a low success rate of computer-aided drug design. The limitations 
of the docking score are obvious. In statistical physics, the free energy is calculated from the partition 
function, which is based on a structural ensemble of numerous structures at a particular temperature. On the 
other hand, the conventional docking score is calculated from a single protein-compound complex structure. 

Many reports have used molecular dynamics (MD) simulations to analyze protein-compound 
docking and to calculate the AG. Even if the protein-ligand complex structure is unknown, ab-initio 
MD docking simulations can be used to reveal the protein-ligand complex structures and the free 
energy landscapes [11-14]. In an explicit water model, if the protein-ligand complex structure is 
known, the binding free energy and the potential of mean force (PMF) along the dissociation path can 
be obtained using the filling potential (FP) method [15], the meta dynamics method [16,17], the 
MP-CAFEE method [18], the smooth-reaction path generation method [19] and Jarzynski's method [20]. 

There have been several reports on the calculation of protein-ligand binding free energy by 
semi-empirical methods, since the ab-initio free energy calculation is still very time-consuming. The 
molecular- mechanical Poisson-Boltzman surface-area (MM-PBSA) method [21], the linear interaction 
energy (LIE) method [22] and the COMBINE method [23~29] have succeeded in reproducing the 
trends of AGs for single target proteins. These methods have been successful in practical use, but the 
parameters used in these methods must be optimized for each target protein. 

We previously proposed a direct interaction approximation (DIA) method for the AG estimation [30]. 
This method estimates the AG value based on the direct interaction between the protein and the ligand 
just as in the COMBINE method, but the weighted parameters for residues are set to fixed values as in 
the LIE method. In the current study, we modified the DIA method in order to use target-independent 
parameters. Since previous authors have introduced a ligand-entropy term in their docking studies [5,6], we 
also examined the ligand-entropy term. In addition, because the mobility of solvent water molecules has 
been analyzed in previous reports [31,32], and we also examined the effect of the solvent water 
mobility herein, but used a different method of analysis for this purpose. 

2. Results and Discussion 

The brief explanation of the previously developed direct interaction approximation (DIA) [30] is 
introduced in Section 2.1. The ligand-entropy term is the first additional term to the original DIA and it 
is introduced in Section 2.2. The stability of hydration shell of each residue is the second additional 
term to the original DIA and it is introduced in Section 2.3. The ligand-entropy term and the stability 
of hydration shell were examined by using the protein-ligand complex structures in Section 2.4. These 
additional two terms improved the accuracy and the physical consistency of the DIA model. These 
results showed that the trajectory average of the protein-ligand interaction improved the estimation of 
the protein-ligand binding free energy. In Section 2.5, we showed the trajectory average of the docking 
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score can improve the binding free energy estimation and the consensus score of the DIA result and the 
docking score improved the correlation between the experimental and the calculated binding free energies. 

2.1. Original Direct Interaction Approximation (DIA) Method 

In our previous study, we developed the DIA method to estimate AG [30]. The fluctuation of the 
accessible surface area (ASA) or the dihedral angles of the system was introduced as the entropy term 
of the AG value, and the estimation accuracy reached 1.5 kcal/mol for several tens protein-compound 
complex structures. Here, we will explain the DIA method briefly. In the original DIA method without 
a direct solvent effect (DIAV), the AG value is estimated as follows [30]: 

AG DIAV = a2_ l <E (i)>xe +P2_ I <E (i)>xe +txS x (1) 

i i 

where E vdW (i) and E ele (i) are the vdW and electrostatic interactions between the i-th residue of the 
protein and the ligand, respectively. Svdw(i) and Sele(i) are the fluctuation of the E vdW (i) and E ele (i), 
respectively. The t*S x term represents the entropy of the system. S x is the fluctuation of a property x. 
In the current study, S x is the fluctuation of the accessible surface area (x = ASA) of the protein- ligand 
complex structure or the all dihedral angles (x = DIH) of the protein over the trajectory. There are five 
parameters: a, a2, [3, [32, and x. 

To represent the van der Waals interaction and the hydrophobic interaction, a Lennard- Jones (LJ) 
8-4-type function has been used instead of the LJ12-6 type function: 

^8-4 = 44 A 8 - A 4 ] , well depth= 8, R e = \[2R 0 (2) 
R R 

where R e is the equilibrium distance. The R e and the well depth values are set to the same values 
obtained from AMBER param99 [33] and the general AMBER force field (GAFF) [34]. The data 
sampling MD simulation is performed with the conventional AMBER force field (LJ 12-6 potential), 
and the DIA analysis is performed with Equation (2). 

In the ligand-binding pocket, the effective dielectric constant (8 e ff) should be different at each point, 
since the s e ff values of proteins are 2-4 and the s e ff of water is 78.5. The E eIe (i) should be scaled by the 
s e ff. Therefore, the electrostatic interactions in the DIAV method were modified and we named the 
modified method as the direct interaction approximation with solvent (DIAS) method [30]: 

AG DMS = aZ < E vdw (i) > x + < E m ;%) > x^* 0 + r x S x (3) 

where E mo d e e (i) is the E e e (i) value scaled by s e jf. The s e ff value could be calculated from the ratio 
between the electrostatic force calculated in the explicit water model and that in a vacuum as follows: 

E eIe (i) 

J b eff 

where Ej ele (i) is the electrostatic interaction between the i-th residue and the j-th atom of the ligand in a 
vacuum. Here, the effective dielectric constant is given by: 
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\(p real .p real } 

M ^ (5) 

c J \ \ ( J7 vac I7 vac \ 
£ eff \ { l , ■ 1 . ) 

where Fj real and Fj Vac are the electrostatic force acting on the j-th atom in the explicit water model and 
in a vacuum, respectively. The F real and F vac were calculated by the molecular dynamics simulation in 
the explicit water model and in a vacuum, respectively. 

The scale factor \le! € could be an unrealisticalry large value when the denominator of Equation (5) is nearly 
zero. Thus, we introduce a parameter x and the scale function as follows: 

s ] eff = 1 + exp( -xxsj ff ) (6) 

This parameter, \ld eS , was used as the scale factor, and the previous study showed that the optimal 
value was 0.6 [30]. Note that the actual effective dielectric constant corresponds to el/P . 

2.2. Intra-Molecular Ligand-Entropy Term 

In the current study, the entropy change of the ligand was taken into account in the AG estimation. 
The rotatable bonds of the ligand can freely rotate in its unbound form, and these bonds can be fixed 
into a single conformation in its bound form. Thus, the entropy of the ligand decreases during the 
protein-ligand binding process. We added the ligand-entropy term (TASii gan d) as follows [5,6]: 

TAS ligand = wxk B TxN rot x]n3 (7) 

where N rot is the number of rotatable bonds (single bonds between heavy atoms) in the non-ring part. 
The number of possible conformers is 3 Nrot without the ring part. Considering the intra-molecule atomic 
collision, the number of conformers can be less than 3 Nrot , and so an additional parameter w is introduced. 

First, the ligand entropy of the ring parts was examined. The number of conformers of a ring part 
was approximated by 2 (Nrot - ring ~ 3) or 3( Nrot - rin § " 3 >, where N rot . ring is the number of rotatable bonds (single 
bonds between heavy atoms) in the ring parts, since the three-membered ring has only one conformer 
and the torsion angles of ring parts are restricted compared to the free rotation. We examined the 
importance of the ring-entropy term by multiple linear regression analysis of the data of 34 protein-ligand 
complexes. Consequently, the ring-entropy term did not improve the estimation accuracy of AG. Thus, 
in the current work, Equation (7) was used as the ligand-entropy term. 

2.3. Hydration Effect of Each Residue of the Target Protein 

In computer-aided drug design, crystal water molecules are often replaced by ligand atoms to design a 
high-activity compound [35]. This is an empirical procedure known among medicinal chemists. The 
amino-acid residues around the crystal water molecules should be important to the protein-ligand 
binding interaction (hot spot). To detect the hot spot, the mobility of water molecules is observed by 
MD simulations. In the current study, the MD simulation of apo protein in water was performed at 
room temperature, and the mobility of the water molecules was observed around the i-th residue. 

The mobility of water is measured by the ligand exchange rate. In the current study, a residue-based 
ligand exchange rate for the i-th residue (<Hj>) was introduced: 
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< H i >=< H™ X — ) / ^ (8) 

m=l 

Here An," 1 is the number of water molecules exchanged at the m-th step of the sampled MD trajectory 
around the j-th atom. The j-th atom belongs to the i-th residue and Ni is the total number of atoms 
(including H atoms) of the i-th residue. Nstep is the total number of the sampled MD steps. The water 

o 

molecule, whose distance to the j-th residue is less than 6 A, is taken into account as the ligand of the 
j-th atom in Equation (8). Since the weight for each energy term of Equations (1) and (3) 
[exp(-a2xS V dw(i)) and exp(-(32xS e ie(i))] corresponds to the probability, the weight for the i-th residue 
is a dimensionless parameter. We assume that the weight of the amino-acid residues is a function of 
<H;>, since it could be a measure of the stability of hydration shell around the i-th residue. The higher 
the value of <H;> is, the more important is the i-th residue is in the protein-ligand interaction. Thus, the 
weight of residues should be a monotonically increasing function of <Hj>. We apply exp(y<H;>) as a 
simple function for approximating the weight, where y is a positive parameter. In the current study, the 
trajectory was sampled every 2 psec and the minimum, maximum and average values of <H;> were 
0.0, 0.16, and 0.042, respectively. These values correspond to 1, 0.26, and 0.96 of the exp(y<H;>) 
values. The average <H;> values in the ligand binding pocket and on the protein surface were 0.0478 
and 0.0488, respectively. The ratios of the residue with <Hj> less than the average <H;> value were 
54% and 48% in the ligand binding pocket and on the protein surface, respectively. 

Water molecules in the bottom of the pocket hardly move and the contact number of the bottom of 
the pocket should be large. Thus there is a correlation between the water mobility and the contact 
number. The contact number is the number of atoms (protein atoms only, excluding the solvent and 

o 

ligand atoms), whose distances from the atom in question are less than 6 A. In this assumption, the 
<H;> value for the i-th residue is estimated instead of Equation (8) as follows: 

Y(C m -C ) 

< H i >=< Hf >= £ (^ ) / Nstep (9) 

m=l N j 

Here Cj m is the contact number of the j-th atom at the m-th step of the sampled MD trajectory and 
C avg is the average value of Cj m . The j-th atom belongs to the i-th residue and Ni is the total number of 
atoms of the i-th residue. In the current study, the correlation coefficient between the <H™> and <H- > was 
0.32. The average (C avg ), minimum and maximum <H- > values were 73.51, 0, and 106.2, respectively. 

2.4. AG Estimation by the DIA Method 

In the current study, we used the following 6 models to examine the ligand-entropy term and the 
hydration effect of residues. 

Model 1: The original DIAV model described in Equation (1). Here, a2 and [32 were set to zero. 
Model 2: The original DIAS model described in Equation (3). 

Model 3: The DIAV_L model, where the ligand-entropy term is added to the original DIAV model: 



AG DIAV L = a^< E vdw (i) > + PYs< > + txS x + TAS 



ligand (JO) 
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Model 4: The DIAV_W model, where the weight of each residue is calculated from the hydration 
solvent water. Here <Hi> = <Hi W > as in Equation (8): 

AG DIAVW = a^expiy < H t » < E vdW (i) > + /?£exp( r < H i » < E e,e (i) > + rxS x (n) 

Model 5: The DIAV_LW model. Here, the ligand-entropy term and the weight for each residue are 
added to the original DIAV model. The weight for each residue is calculated from the hydration solvent 
water. Here <H> = <H W > as in Equation (8): 

^G DIAVLW = a^pir < H i » < E vdW {i) > + /?Xexp( r < H. » < E ele (i) > + rxS x + TAS [igand {n) 

Model 6: The DIAV_LC model. The weight for each residue is estimated based on the contact 
number. The model equation is Equation (12) with the relation <Hi> = < H- > as in Equation (9). 

Table 1 shows the computational average error and the correlation coefficient between the 
experimental values and the values calculated by these six models. The results were obtained by 
leave-one-out cross-validation tests. In the leave-one-out cross-validation test, one data is selected as 
the test data that is to be predicted and the other data are used as the teaching data to generate the 
prediction model equation. The test data is selected one after another in the given data set until all data 
are selected as the test data. The property x of Sx (entropy term) was fixed to x = DIH (the fluctuation 
of the dihedral angles), since the AG accuracy obtained by x = DIH was better than that obtained by 
x = ASA (the fluctuation of the accessible surface area). The y values were optimized to minimize the 
AG estimation error, and these values were set to -6.115 and 0.00613 for the DIAV_LW and 
DIAV_LC methods, respectively. 

Comparing the average error obtained by the DIAV and DIAV_L models, the consideration of the 
ligand-entropy term improved the accuracy. In addition, comparing the average error obtained by the 
DIAV and DIAV_W models, the consideration of the weight of the residues improved the accuracy. 
The combination of both the ligand-entropy term and the weight of the residues improved the accuracy 
(DIAV_LW). Among the six models examined, the DIAV_LW model showed the best accuracy and 
the DIAV_C model showed the second best accuracy. The accuracy of the DIAV_L model was almost 
the same as that of the DIAS model. Since the number of parameters of the DIAV_L model was 
smaller than that of the DIAS model, the second best model should be the DIAV_L model. 

Figure 1 shows the correlation between the experimental and calculated AG values obtained by the 
DIAV, DIAV_LW and DIAV_LC methods. These values were obtained by the molecular dynamics 
simulation starting from the experimentally determined protein-ligand complex structures. 
The computational detail was described in the data preparation section. It is clear that the 
DIAV_LW/DIAV_LC methods gave a better correlation than the DIAV method. 
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a: standard deviation of the difference between calculated and measured binding free energy, b: We also 
applied the multiple linear regression (MLR) to the 34 protein -ligand complex data. "Average Error (MLR)" 
is the averaged error obtained by the MLR. The error of the MLR is always smaller than the error obtained by 
the cross-validation test. 
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Figure 1. Cross-validation results obtained by the DIAV, DIAV_LW and DIAV_LC methods: 
The experimental data (AG exp ti) and the calculated values (AG ca i c ). Open circles, filled circles 
and filled triangles represent the results obtained by the DIAV, DIAV_LW and DIAV_LC 
methods, respectively. 
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Table 2 shows the average, deviation, and minimum and maximum values of the optimized 
parameters (a, P, x and w) of Models 1-6 of the 34 cross-validation tests. The % of the negative values is 
also summarized. The other parameters, i.e., a2, [32, and y, are omitted. The smaller the deviation of the 
parameter is, the less dependent on the target protein the model is. In particular, the sign of parameter [3 is 
important. Negative values of [3 are physically unreasonable. Namely, negative P implies that repulsive 
ligand-protein interactions stabilize the free energy of binding. In the DIAS and the DIAV models, the 
P value was negative in 2.86% of the cases (one model among total 34 cross-validation test models). In 
contrast, the average a, x and w values were almost identical among the DIAV_L, DIAV_LW, DIAV_LC 
and DIAS models. The consideration of the ligand-entropy term (DIAV_L model) slightly improved 
the problem of the negative P parameter. In addition, the weight of residues (DIAV_W model) also 
slightly improved the problem of the negative P parameter. In the DIAV_LW and DIAV_LC models, 
all P parameters were positive in the 34 cross-validation tests. 

Among these six models, the deviation of the DIAV_LW model was the smallest. Considering the 
average error and the deviation of the parameters, the DIAV_LW model was the best of the 6 models 
and the DIAV_LC model is the second best. In the drug design, the prediction accuracy of unknown 
compounds for a new target protein is more important than the regression of the activities of known 
active compounds for a known target. Thus, the parameters of the computational model must not 
depend on the target proteins. From this point of view, the DIAV_LW or DIAV_LC model is desirable. 



o 
M 



< 



-2 
-4 
-6 
-8 
■10 
■12 
■14 
-16 



Pharmaceuticals 2013, 6 612 



Table 2. Summary of parameters determined by the cross-validation tests. 



DIAV 


a 


P 


X 


w 


Average 


a r\i a 1 n 1 n 
U.U341 / 19 


A AA1 

U.UU1 033 


A AAAO 1 AO 


A AAAAAAA 

U.UUUUUUU 


Deviation (a) 


A AAAC4AC 


A AA1 1 OH A 

U.UU1 lo /4 


A AAAAA07 

U.UUUUUo / 


A AAAAAAA 

U.UUUUUUU 


Min 


A A^ C1 1 

U.U323M 1 


A AAO O O AT 

— U.UU3ooU/ 


A AAAO A 1 O 

— U.UUU243© 


A AAAAAAA 

U.UUUUUUU 


Max 


A A^3 C7C£/| 

U.U3:> Oo4 


a AA/imno 
U.UU49/VO 


A AAAOAOT 

— vj.vjvjvjIUZi 


A AAAAAAA 

U.UUUUUUU 


Negative value 


0.0000000 


0.0285714 


1.0000000 


0.0000000 


1~\T A 1 7 T 

D1AV_L 


a 


0 

P 


T 


w 


Average 


A A/^ 

U.U3 /UlVo 


A AAOA^C 1 

U.UU296M 


A AAAAACA 

— U.UUUUUdU 


A 1 T/IA1 £.C\ 

U. 1 /4VloV 


Deviation (a) 


A aaat<oo 


A AAAO A <C\ 

U.UUUo4:>U 


A AAAAAAO 

(j.vvvvvvl 


A A1 AACn 1 

U.U19UD21 


Min 


U.U3DUV33 


A AAAA/^ A 1 

— U.UUUUo41 


A AAAAAC4 

— U.UUUUUM 


A 1 1 QOQOQ 

U. 1132383 


Max 


U.U3Vo24y 


U.UU4/2U4 


A AAAAA/1 C 

— U.UUUUU4D 


A O CAT /I 

U.232DW4 


Negative value 


rv (\(\(\(\(\(\(\ 

0.0000000 


0.0285714 


1.0000000 


0.0000000 


D1AV_W 


a 


P 


T 


w 


Average 


U.U34oo23 


a aao 1 non 
U.UU21V2V 


A AAAOAC-1 

— U.UUU2UM 


A AAAAAAA 

U.UUUUUUU 


Deviation (a) 


A AAAC3 O O 

U.UUU:>3oo 


A AA1 1 C\1 £. 

U.UU1 1U36 


A AAAAAOQ 

U.UUUUU83 


A AAAAAAA 

U.UUUUUUU 


Min 


U.U3292/3 


A AAOAO /I O 

— U.UU3U242 


A AAAIOAA 

—v.vjvjvjllyv) 


A AAAAAAA 

U.UUUUUUU 


Max 


U.U3o2o9V 


A AACAAAC 


A AAA1 070 

— U.UUUlo / o 


A AAAAAAA 

U.UUUUUUU 


Negative value 


0.0000000 


0.0285714 


1.0000000 


0.0000000 


D1AV_LW 


a 


0 

P 


T 


w 


Average 


A A/1 1^1 £.1 

U.U413163 


A AA^IAOO 

U.UU62U33 


A AAAAA/^T 

— U. UUUUUo/ 


U.153ollo 


Deviation (a) 


A AAAT2 OO 


A AAA7AA7 
U.UUU/9U/ 


A AAAAAAA 

U. UUUUUU2 


A A1 A C\C\A A 
U.U 14UU4U 


Min 


U.U3V20/ / 


A AAO /: 

U.UU342 16 


A AAAAA7 1 

— U. UUUUU / 1 


A 11C/IA/I/I 

U. 1254U44 


Max 


U.U43448U 


A AAC7 1 /^O 

U.UUo /lo2 


A AAAAA/iQ 

— U. UUUUU63 


A 1 C\A A A AH 
U. 194444/ 


Negative value 


A AAAAAAA 

0.0000000 


A AAAAAAA 

U.UUUUUUU 


1 AAAAAAA 

1.0000000 


A AAAAAAA 

0.0000000 


L)1AV_LC 


a 


0 

p 


T 


w 


Average 


A A^2 /1 1 A/I 

0.034304b 


U.UU42VD0 


A AAAAAAA 

—o.ooooo/o 


0. 1 143zyj 


Deviation (a) 


A AAA£ 1 OA 

o.ooooizy 


A AA1 1 ACA 

U.UU1 IVjU 


A AAAAAAO 

0.0000002 


A AAmo 1 £ 

o.ooy /zio 


Min 


A A^O 1 ^70 

0.03213 /8 


A AAATOQ C 

U.UUU2o3d 


A AAAAAT^ 

—0.00000/0 


A AO ylOOOC 

0.0942835 


Max 


A A^2/£^2AA1 

O.OJoJOOl 


U.UUVUDdd 


A AAAAA/£^7 

— O.OOOOOo/ 


A 1 A 1 /I <CA/1 

U. 14145U4 


Negative value 


A AAAAAAA 

u.uuuuuuu 


A AAAAAAA 
U. UUUUUUU 


1 AAAAAAA 

1. UUUUUUU 


A AAAAAAA 

U.UUUUUUU 


DIAS 


a 


P 


T 


w 


Average 


0.0392333 


0.0030804 


-0.0000053 


0.0000000 


Deviation (a) 


0.0005573 


0.0010426 


0.0000002 


0.0000000 


Min 


0.0375654 


-0.0017236 


-0.0000056 


0.0000000 


Max 


0.0409116 


0.0055269 


-0.0000049 


0.0000000 


Negative value 


0.0000000 


0.0285714 


1.0000000 


0.0000000 



The consideration of the ligand entropy and the weight of each residue did not improve the DIAS 
model. In the DIAS model, the weight of each residue is already considered by using the parameters 
a2 and [32. Thus, the newly introduced weight with the y parameter would count the weight of the 
residue twice. Considering the number of parameters (a, |3, x, y and w in the DIAV_LV/DIAV_LC models, 
and a, [3, x, a2, [32 and y in the DIAS model), the DIAV_LW/DIAV_LC models have a smaller number of 
parameters than the DIAS model. Since a model with a small number of parameters should, in 
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principle, be better than that with a large number of parameters, the DIAV_LW/DIAV_LC models are 
better than the DIAS model. 

The trajectory dependence of the models was examined. The above results were obtained from the 
2-nsec trajectories. Figure 2 shows the time dependence of the DIAV, DIAV_LW and DIAV_LC 
results. When the 1-nsec trajectories were used, the results were slightly worse than those shown 
above, but the difference was not statistically significant. In the DIAV model, the P value was 
negative in 5.71% of the cases (lett and lhsl). In the DIAS, DIAV_L and the DIAV_W models, the P 
value was negative in 2.86% of the cases (lett in the all models) and no negative P value was observed 
in the DIAV_LW/DIAV_LC models, just as in the above results. The lett structure is thrombin, but 
the other thrombin structures (letr, lets) did not show the problem. Currently, the reason of the 
problem is unclear. 

Figure 2. Time dependence of the DIAV, DIAV_LW and DIAV_LC results. Filled circles, 
open squares and open circles represent the results obtained by the DIAV, DIAV_LW and 
DIAV_LC methods, respectively. The averaged error is plotted vs. the sampling-time length. 




We also estimated the binding free energies of non-active compounds, since evaluation of the binding 
free energies of the both active and non-active compounds are essential in practical use. 
We docked three GPCR (G-protein coupled receptor: membrane protein) ligands to the proteins and 
calculated the binding free energies of them by the DIAV_L method. These compounds were alprenolol 
(a P-adrenergic inverse agonist), fenoterol (a P-adrenergic inverse agonist) and cetirizine (Hi receptor 
inverse agonist). They are non- active compounds of the proteins, because there are no GPCR in the 
target proteins. The condition of the MD simulation was the same as those used in Table 1. The 
parameters of the DIAV_L method were determined for the MLR of the all 34 proteins used in Table 1. 
The results are summarized in Table 3. In some cases, the pocket sizes of the ligand-binding sites of the 
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proteins were too small for the GPCR ligands to be docked, and so MD simulations were not possible 
for them. Obviously, such ligands should be non-active compounds. The average values of the 
estimated binding free energies for the non-active compounds were weaker than those of the original 
ligands. In all 19 cases, the AGs of alprenolol were weaker than those of the original ligand. The AGs 
of fenoterol and cetirizine were weaker than those of the original ligand in 69% and 56% cases, 
respectively. For some proteins, feneterol and cetirizine show stronger affinities estimated than those 
for the original ligands, but most of their absolute bindings are weak. 

Table 3. The binding free energies estimated for the non- active compounds (kcal/mol). 



PDB ID 


AG ,i 


Original 
i^iiiiii 

li^uiid 


Alprenolol 


Difference a 


Fenoterol 


Difference a 


Cetirizine 


Difference a 




— Q 57 


O . V./V' 


-6 85 


-1.22 


-8.21 


0 14 


-9 28 


1.22 


1 ahfl 


-7 ^0 


-8 40 


-6 13 


-? 97 


-6 7? 


-1 


/ .yj 


-0 47 

\J.'-T 1 


1 nnn 


-10 50 


-1 1 63 


-2 77 


-8 86 


-4 50 


-1 12 


-5 69 


-5 93 

u .y ~j 


1 chv 




-8 8Q 


— S 84 


04 


-7 51 


-1 ^8 
1 .JO 


O.Ju 


-0 58 

W..JO 


ldog 


-5.48 


-9.05 


-5.18 


-3.87 


-7.75 


-1.30 


-5.08 


-3.97 


ldwb 


-3.98 


-5.45 


-5.44 


-0.01 


-6.56 


1.11 


-8.24 


2.80 


lebg 


-14.76 


-6.74 


0.00 


-6.74 


0.00 


-6.74 


0.00 


-6.74 


lepo 


-10.85 


-14.42 


-5.64 


-8.78 


-7.30 


-7.12 


-8.49 


-5.93 


lrbp 


-9.17 


-8.76 


N.D. b 


N.D. b 


N.D. b 


N.D. b 


-8.69 


-0.08 


lstp 


-18.27 


-6.59 


N.D. b 


N.D. b 


N.D. b 


N.D. b 


-5.96 


-0.63 


ltnh 


-4.59 


-5.59 


-4.39 


-1.20 


-5.62 


0.03 


-6.13 


0.54 


lulb 


-7.23 


-6.19 


-5.45 


-0.74 


-6.23 


0.04 


-8.98 


2.79 


2gbp 


-10.36 


-10.14 


-7.16 


-2.98 


-8.74 


-1.40 


-10.24 


0.11 


2ifb 


-7.41 


-8.60 


-5.81 


-2.79 


-7.09 


-1.51 


-9.01 


0.41 


2tsc 


-11.62 


-8.23 


-5.68 


-2.55 


-6.48 


-1.75 


-8.69 


0.47 


2ypi 


-6.58 


-6.92 


-4.68 


-2.24 


N.D. b 


N.D. b 


N.D. b 


N.D. b 


3ptb 


-6.46 


-4.96 


-4.49 


-0.48 


-5.89 


0.93 


-5.64 


0.68 


4dfr 


-13.22 


-8.42 


-5.16 


-3.26 


-5.64 


-2.79 


-6.66 


-1.76 


6cpa 


-15.71 


-11.68 


-6.82 


-4.86 


-7.77 


-3.91 


-9.75 


-1.93 


Average 


-9.57 


-8.35 


-5.15 


-3.29 


-6.38 


-2.15 


-7.38 


-1.06 



a: the energy difference between the calculated AG of the original ligand and the non-active ligand.b: Not 
Determined, because the pocket sizes of the ligand-binding sites of the proteins were too small for these 
ligands to be docked, and so MD simulations were not possible for them. 



The whole protein set included four thrombins, three HIV-1 proteases and five trypsins. We have 
examined the Spearman's rank correlations for these ligands of the same target proteins. The parameters of 
the DIAV_L, DIAV_LW and DIAV_LC methods were determined based on the 22 proteins excluding 
these 12 target proteins. The AGs, error of the AGs, and the correlation coefficients are summarized in 
Table 4. The total number of HIV-1 proteases was 11, since we added eight new data points. These results 
obtained by the DIAV_L, DIAV_LW and DIAV_LC methods are similar to each other. The trends of 
the AGs are almost correctly reproduced. 
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Table 4. AG values (kcal/mol) of the same target proteins and Spearman's rank correlations. 



Thrombin 


AG expt i 


AGdiav l 


Error 


^Gdiav lw 


Error 


^Gdiav lc 


Error 


ldwb 


-3.98 


-5.15 


1.17 


-5.02 


1.04 


-5.57 


1.59 


lett 


-8.44 


-9.9 


1.46 


-9.74 


1.31 


-9.81 


1.37 


letr 


1 A AA 

-10.09 


A A 

-9.9 


0.19 


A OA 

-9.89 


A O 

0.2 


1 A T"> 

-10.22 


A 1 1 

0.13 






-10 Q 
Lvj.y 


0.72 


-10 7fi 


0.86 


-10 46 


1.16 


Averaged error 
(kcal/mol) 


- 


- 


0.89 


- 


0.85 


- 


1.06 








l.Ul 




V.yJ 




1 .zu 


Correlation 
coefficient 


- 


- 


0.97 


- 


0.97 


- 


0.96 


Spearman s rank 


- 




1 




1 


- 


1 


correlation 




- 


- 




H1V-1 Protease 


A /~* 


A*JDIAV_L 


Error 


A*JDIAV_LW 


Error 


A*J£)IAV_LC 


Error 


lk6p 


-8.84 


-11.71 


2.87 


-11.74 


2.90 


-11.78 


2.94 


lajv 


-10.59 


-10.36 


0.23 


-10.39 


0.20 


-10.13 


0.46 


lajx 


-10.86 


-9.89 


0.97 


-9.91 


0.95 


-9.68 


1.18 


lhih 


-10.97 


-11.67 


0.70 


-11.67 


0.70 


-11.73 


0.76 


lhtf 


-11.04 


-11.57 


0.53 


-11.59 


0.55 


-11.86 


0.82 


laaq 


-11.45 


-13.15 


1.70 


-13.13 


1.68 


-12.96 


1.51 


lhpv 


-12.57 


-12.79 


0.22 


-12.87 


0.30 


-13.06 


0.49 


lhvr 


1 O f\'~1 

-12.97 


-14.79 


1 O O 

1.82 


1 A AO 

-14.93 


1 A/' 

1.96 


1 A f C 

-14.65 


1.68 


lhvk 


-13.79 


-13.63 


0.16 


-13.65 


0.14 


-13.70 


0.09 


lv J 


1 A 

-14. oz 


1 1 oo 

-12.82 


1 OA 

1.80 


-12.85 


1 T7 

Y.I 1 


1 o on 
-12.89 


1 HI 

Y.I 5 


ldif 


-14.63 


-13.76 


0.87 


-13.77 


0.86 


-13.82 


0.81 


Averaged error 
(kcal/mol) 


- 


- 


1.08 


- 


1.09 


- 


1.13 


oJJ 






1. JO 




1 T7 
1. J / 




1 in 
1.3 / 


Correlation 
coefficient 


- 


- 


0.68 


- 


0.67 


- 


0.68 


Spearman's rank 
correlation 


- 


- 


0.78 


- 


0.75 


- 


0.81 


Trypsin 


AG expt l 


4Gdiav_l 


Error 




Error 


4Gdiav_lc 


Error 


ltng 


-4.00 


-5.45 


1.45 


-5.36 


1.37 


-2.69 


1.31 


ltnn 


A CC\ 

-4.59 


c o a 

-5.29 


0.70 


-5.20 


0.61 


-5.50 


A A 1 

0.91 


3ptb 


-6.46 


-4.92 


1.54 


-4.83 


1.63 


-5.15 


1.31 


lpph 


O A O 

-0.48 


O 1 o 

-8.32 


A 1 £ 

O.lo 


O ^ A 

-8.30 


A 1 O 

0.18 


O C 1 

-8.51 


A AO 

0.02 


lppc 


-8.80 


-9.32 


0.52 


-9.31 


0.51 


-9.53 


0.72 


Averaged error 
(kcal/mol) 






0.88 




0.86 




0.86 


SD a 






1.03 




1.02 




0.98 


Correlation 
coefficient 






0.86 




0.86 




0.93 


Spearman's rank 
correlation 






0.60 




0.60 




0.90 



a: standard deviation of the difference between calculated and measured binding free energy. 
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2.5. Consensus Score with the Trajectory Average of the Docking Score 

Next, we examined the trajectory average of the docking score. The Sievgene protein-compound 
docking program was used to calculate the protein-ligand docking score [7]. The trajectory average 
improved the correlation between the experimental binding free energy and the averaged docking 
score. Namely, the correlation coefficients between the experimental binding free energy and the 
averaged docking score were 0.751 and 0.745, with and without the trajectory average, respectively. 
The actual docking scores calculated by three different programs were summarized in the supplementary 
data. On the contrary, the DIAV_L (R = 0.76) and DIAV_LW (R = 0.78) methods gave R = 0.76 and 0.78, 
respectively, slightly better than those for the averaged docking scores. However, the differences between 
their estimated binding free energies and the experimental ones were 1.39 kcal/mol and 1.33 kcal/mol, 
respectively. They are much smaller than those by the averaged docking scores, 1.63 kcal/mol and 
1.89 kcal/mol with and without trajectory average, respectively. We must note that the results should 
strongly depend on the ensemble generated by the MD simulation [36]. Our method should be applied 
the correct protein-ligand complex structures or the protein-ligand complex structures must be the 
equilibrium states, otherwise the calculated AG values should drift depending on the simulation time. 

The consensus score of the DIA model and the docking score was also examined. The simple sum 
of the AG value obtained by the DIAV_LW and Sievgene docking score gave a correlation coefficient 
between the consensus score and the experimental AG of 0.796. The simple sum of the AG value 
obtained by the DIAV_LC and Sievgene docking score gave a correlation coefficient between the 
consensus score and the experimental AG of 0.782. Thus, the consensus score worked well and it slightly 
improved the AG estimation. 

Recently, a machine-learning approach was applied to improve the docking score. Such new method 
showed the AG standard deviation (SD) error of 1.5 kcal/mol and the correlation coefficients between 
the experimental binding free energy and the docking score reached 0.76 based on a single structure [37]. 
This result is better than our current result (R = 0.75-0.76, SD = 1.7-1.8 kcal/mol; see Table 1), however 
the used protein-ligand datasets were different to each other. The accuracy of the other docking score 
could be improved considering the suitable ensemble of the protein-ligand complex structures. 

3. Method: The Docking Score Calculation 

A protein-compound docking simulation was performed by the program Sievgene, which is a 
protein-ligand flexible docking program for in silico drug screening [7]. This program generates many 
conformers (100 conformers by default) for each compound and keeps the target protein structure 
rigid, but with soft interaction forces adapting its slight structural change to some extent. The Sievgene 
scoring function was designed to consider the structural change of the target protein. In the inner 
region of the target protein, the protein is approximated as an elastic body, while the atomic pair-wise 
scoring function is applied in the outer region of the target protein. This docking program was developed 

o 

with a performance yielding about 50% of the reconstructed complexes at a distance of less than 2 A 
RMSD for the 132 complexed receptors with the compounds in PDB. The results predicted by our 
program were almost the same as those predicted by other docking programs [7]. The docking score 
(Hdock) to estimate the protein-ligand binding free energy was determined as: 
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^dock C rot ' ^ rot C AV ' ^ ASA ^vdw) C ele ' ^ele C hyd ' ^hyd C mtra-vdW ' ^"mtra-vdW (15) 

where N rot , E A sa, E vd w, E e i e , E hyd , and Ei ntra -vdw represent the number of rotatable bonds of the ligand 
molecule, the hydrophobic energy due to the accessible surface area, the vdW energy, the protein- ligand 
Coulombic potential, the hydrogen bond energy, and the intra-molecular vdW energy of the ligand for 
Sievgene [7]. Also, c rot , c A v, c e u, Ch y d, and Cmtm-vdw are the optimized coefficients for each energy term. 
For each atom type, the sum of Easa and E v dw gives a grid potential, and both energy terms are always 
simultaneously calculated. Thus, these two terms share the same coefficient, cav- Sievgene utilizes the 
grid potential to calculate each energy term except for the intra-molecular interactions. In this study, a 
mesh size of 60 x 60 x 60 was adopted. 

4. Data Preparation 

To determine the coefficients for the AG scores for several current models, we performed a 
protein- ligand docking simulation based on the known complex structures registered in the Protein 
Data Bank. The data and the procedure were almost the same as those used in the previous study [30]. 
Here, 34 complexes accompanied by the experimental binding free-energy values were selected from 
the database that was used to determine the AG scores of the PRO_LEADS [6]. The PDB identifiers, 
the names (protein names and ligand names), the molecular weights (MW), the number of hydrogen 
bond acceptors (HA) and the number of hydrogen bond donors (HD) of ligands are summarized in 
Table 5. There was no peptide- like compound. The MWs were distributed from 114 to 606 Da. To assess 
the ligand diversity, we calculated the average Tanimoto index and the standard deviation of the all 
34 ligands x all 34 ligands by using Maximum Common Substructure (MCS) method [38]. The average 
Tanimoto index and the standard deviation were 0.29 and 0.19, respectively. These values showed that the 
used ligand molecules were diverse. In the test dataset, the metalloproteins were removed from the 
present analysis. Metal atoms (Zn and Fe atoms) formed covalent bonds with O and S atoms of the 
ligands, and the classical force field that we applied could not represent the covalent bond. Thus, the 
present method cannot calculate AG values for metalloproteins with high precision. 

The structural ensembles generated from the PDB structure given by MD in explicit water were 
prepared as follows. All target proteins were prepared with ligands (forming a protein-ligand complex 
structure). In the previous study, all metal atoms in the systems were removed, since the target proteins 
were not metalloproteases. Some non-metalloproteins include metal atoms those bind to the proteins. 
In the current study, all metal atoms of the PDB files were included in the MD simulations. The force 
fields and the charges of the protein atoms originated from AMBER parm99 [36]. The atomic charge 
of each ligand was determined by the restricted electrostatic point charge (RESP) procedure using 
HF/6-31G*-level quantum chemical calculations [39]. We used Gaussian98 to perform the quantum 
chemical calculations [40]. The initial coordinates of protein and ligand molecule of each data were 
fixed to the experimentally determined coordinates. The whole structure of each protein was embedded 
in a sphere of TIP3P [41] water molecules (CAP water), including ion particles of 0.1% Na + and CP, 
in order to neutralize the total charge of the systems. The center of the sphere was set at the mass 
center of the protein. The shortest distance between the protein atom and the CAP sphere wall was 
set to 10 A. 
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Table 5. List of the proteins and ligands used. 



PDB ID 


Protein 


Ligand 


MW 


HA 


HD 


labe 


L-arabinose-binding protein 


L-arabinose 


150.1 


5 


4 


labf 


L-arabinose-binding protein 


D-fucose 


161.2 


5 


4 


lapu 


acid proteinase (penicillopepsin) 


pepstatin 


485.7 


6 


4 


ldbb 


Fab' fragment 


progesterone 


314.5 


2 


0 


ldbj 


Fab' fragment 


aetiocholanolone 


290.4 


2 


1 


ldog 


glucoamylase 


deoxynojirimycin 


163.2 


4 


5 


ldwb 


thrombin 


benzamidine 


120.2 


0 


2 


lepo 


endothia aspartic proteinase 


N-carbonylmorpholine 


131.1 


5 


6 


letr 


thrombin 


MQPA 


509.2 


5 


5 


lets 


thrombin 


NAPAP 


522.6 


4 


4 


lett 


thrombin 


4-tapap 


429.6 


3 


3 


lhpv 


HIV-1 protease 


amprenavir 


505.6 


6 


3 


lhsl 


Histidine-binding protein 


Histidine 


156.2 


3 


2 


lhtf 


HIV-1 protease 


GR 126045 


574.7 


4 


5 


lhvr 


HIV-1 protease 


XK263 


606.8 


3 


2 


lnsd 


neuraminidase 


neuraminic acid 


290.2 


8 


5 


IpgP 


6-phosphogluconate dehydrogenase 


6-phosphogluconic acid 


276.1 


10 


4 


lphg 


cytochrome P450 


metyrapone 


226.3 


3 


0 


lppc 


trypsin 


Napap 


533.6 


4 


4 


lpph 


trypsin 


3-Tapap 


429.6 


3 


3 


lrbp 


retinol-binding protein 


retinol 


286.5 


1 


1 


ltng 


trypsin 


aminomethylcyclohexane 


114.2 


0 


1 


ltnh 


trypsin 


4-fluorobenzylamine 


126.2 


0 


1 


lulb 


purine nucleoside phosphorylase 


guanine 


151.1 


3 


3 


2cgr 


Igg2b (KAPPA) Fab fragment 


guanidineacetic acid 


384.4 


3 


3 


2gbp 


D-galactose / D-glucose-binding 
protein 


D-glucose 


180.2 


6 


5 


2ifb 


intestinal fatty acid binding protein 


palmitic acid 


256.4 


2 


0 


2phh 


P-hydroxybenzoate hydroxylase 


P-hydroxybenzoate 


138.1 


3 


1 


2r04 


rhinovirus 14 (HRV14) 


W71 


342.4 


5 


0 


2tsc 


thymidylate synthase 


1 0-propargyl-5 , 8 -dideazafolic 
acid 


477.5 


7 


3 


2ypi 


triose phosphate isomerase 


2-phosphoglycolate 


156.0 


6 


0 


3ptb 


trypsin 


benzamidine 


120.2 


0 


2 


4dfr 


dihydrofolate reductase 


methotrexate 


454.4 


9 


3 


5abp 


L-arabinose-binding protein 


D-galactose 


180.2 


6 


5 



MW: Molecular weight (Da); HA: Number of hydrogen bond acceptors; HD: Number of hydrogen bond donors. 
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Before an MD calculation was performed for the entire system, an MD calculation for only the 
solvent parts (solvent water and counter ions) was performed with the protein, ligand, and metal ion 
coordinates fixed, so as to bring the solvent parts sufficiently close to an equilibrium state. The 
SHAKE method was used to constrain covalent bonds between heavy and hydrogen atoms in any 
molecule in the system [42]. MD simulations of the entire system were performed using 2.0 fs time 
steps with the temperature set at 310 K; the fast multipole method [43] was used to calculate the 

o 

Coulombic interaction. The cutoff distance of the van der Waals interaction was 12.0 A. The MD 
simulations were performed by using cosgene/my Presto [15]. After equilibration steps of 1,000 ps, the 
protein coordinates were sampled every 2 ps. Finally, we obtained 1,000 structures for each 
target protein in the 2,000 ps production run. The software program myPresto version 4 
(http://presto.protein.osaka-u.ac.jp/myPresto4/index_e.html) was used for the simulation. The 2-nsec 
MD simulations cost average 79 h (max 229 h) using 4-core parallel computation on intel Xeon 5600. 
The trajectory analysis for the DIA method cost average 580 second using single core on intel Xeon 5600. 

5. Conclusions 

We have developed new computational models for protein-ligand binding free energy estimation. 
The DIAV_LC and DIAV_LW models were based on the trajectory average of the protein-ligand 
interaction with the ligand-entropy term. The mobility of the water molecules in the ligand-binding 
pocket was used to calculate the weight of the each residue-ligand interaction of the target protein. The 
interactions of residues around the low-mobility water were weighted comparing to the interactions of 
other residues. The consideration of the ligand entropy and the weight of the residues reduced the 
target-protein dependence of the parameters of the DIA models and consequently the accuracy was 
improved. The average error of AG estimation was 1.3 kcal/mol and the correlation coefficient between 
the experimental values and the calculated values was 0.75, when the correct protein-ligand complex 
structures were provided. The trajectory average of the docking score improved the correlation 
between the docking score and the experimental AG values. In addition, the simple sum of the AG 
value obtained by the DIAV_LW/DIAV_LC and a docking score showed the correlation coefficient 
between the consensus score and the experimental AG of 0.8. Thus, the consensus score worked well, 
and it slightly improved the AG estimation. 
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